% Validation process for NetCDF NEMO model outputs.
%
% --- Call: nemo_validation
%
% --- Input:
%
% --- Output:
%
% Warning: This code is a research code that works in a specific case. This
% is distributed as a shape for future developments on model validation.
% -- Contact for more information: guillaume.charria@ifremer.fr
%
% G. Charria (02/2008)
% ------------------------------------------------------
close all
clear all


% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% ---- READING ...
% ------------------------------------------------------

% ------------------------------------------------------
% -- Reading observations
%
% --- Datasets --- Kind --- Parameter ------------------
%    cruise        woce      T S NO3          
%                  amt       T S CHL +?
%    mooring       bats      T S NO3 CHL PON DON PP
%    climato       woa05     T S NO3
%                  mld       MLD
%    satellite     seawifs   CHL or OWP
%                  eke       EKE
%                  sst       T
%                  pp        PP
%    floats
% ------------------------------------------------------
dataset='climato';
kind='woa05';
parameter='NO3';
cruise=''; % Default
%cruise='amt4'; % Cruise ID
year='1997';
% -- Warning for year: it should correspond to the year of the given
% cruise or mooring


bbio='SBIO3.1';
%bbio='BBIO2.5';
bbio_dir=['/media/disk/NEMO_OUTPUT/',bbio,'/'];

lec=0;
stat=0;
fig=1;

% -----

if (lec==1)

    % dat: lon,lat,time
    dat=r_data(dataset,kind,parameter,cruise,year);
    % (the time is given in MATLAB Serial Number generated using 'datenum.m')

    % ------------------------------------------------------
    % -- Reading model output (depending the observations)
    % ------------------------------------------------------
    % - temporarly test because sst sattelite doesn't exist.
    if (strcmp(kind,'sst'))
        return;
    end
    model=r_model(dat,dataset,parameter,bbio,bbio_dir,year);

    % ------------------------------------------------------
    % -- Interpolation of model and data on a same grid
    % ------------------------------------------------------
    [nmodel,dat]=i_model(dataset,kind,dat,model,parameter);

    % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ End of reading ~~~~


    % ---- Sauvegarde des champs regrilles pour reutilisation plus rapide.
    save([bbio_dir,bbio,'_',dataset,'_',kind,'_',parameter,'_',cruise,'_',year,'.mat'],'nmodel','dat')

else
    
    % ----- Chargement du .mat
    load([bbio_dir,bbio,'_',dataset,'_',kind,'_',parameter,'_',cruise,'_',year,'.mat'])
    
end


% ------------------------------------------------------
% ---- STATISTICS
% ------------------------------------------------------
if (stat==1)
    
    % -- Correction manuelle de points aberrants lies a l'interpolation
    switch kind
        case 'woa05'
            switch parameter
                case 'T'
                    nmodel.data(nmodel.data==0)=NaN;                    
                    nmodel.data(nmodel.data>1e3)=NaN;
                    nmodel.data(nmodel.data<-1e3)=NaN;
                case 'S'
                    nmodel.data(nmodel.data==0)=NaN;                    
                    nmodel.data=nmodel.data(2:end,:,:,:);
                    nmodel.depth=nmodel.depth(2:end);
                    dat.data=dat.data(2:end,:,:,:);
                    dat.depthdata=dat.depthdata(2:end);
                    %nmodel.data(nmodel.data>1e3)=NaN;
                    nmodel.data(nmodel.data<30)=NaN;
                case 'NO3'
                    nmodel.data(nmodel.data==0)=NaN;                    
                    nmodel.data=nmodel.data(2:end,:,:,:);
                    nmodel.depth=nmodel.depth(2:end);
                    dat.data=dat.data(2:end,:,:,:);
                    dat.depthdata=dat.depthdata(2:end);
                    
            end 
        case 'woce'
            switch parameter
                case 'S'
                    switch cruise
                        case 'ar12_g'
                            nmodel.data(nmodel.data<35)=NaN;     
                        case 'a25'
                            nmodel.data(nmodel.data<34)=NaN;     
                            dat.data(dat.data<34)=NaN;     
                        case 'a20'
                            nmodel.data(nmodel.data<30)=NaN;       
                        case {'a22','a02b'}
                            nmodel.data(nmodel.data<32)=NaN;    
                    end
                case 'NO3'
                    switch cruise  
                        case {'a22','a02b'}
                            nmodel.data(nmodel.data==0)=NaN;                       
                    end
            end 
        case 'amt'
            switch parameter
                case 'S'
                    switch cruise
                        case {'amt4','amt5'}
                            nmodel.data(nmodel.data<30.4)=NaN;                    
                    end
                case 'T'
                    switch cruise
                        case {'amt4','amt5'}
                            nmodel.data(nmodel.data==0)=NaN;                    
                            
                    end
            end  
        case 'seawifs'  
            switch parameter
                case {'CHL','OWP'}
                   nmodel.data(nmodel.data==0)=NaN; 
                   %nmodel.data(nmodel.data>1)=NaN;
                   nmodel.data(:,:,1)=NaN;
            end          
            
    end       
    
    dat
    nmodel 
    [sigmodel,sigdata,corr]=s_overview(dat,nmodel);
    %corr=abs(corr);
    disp(['Sigma_model=',num2str(sigmodel),' // Sigma_data=',num2str(sigdata),' // Correlation=',num2str(corr)])
    save([bbio_dir,'Stat_',bbio,'_',dataset,'_',kind,'_',parameter,'_',cruise,'_',year,'.mat'],'sigmodel','sigdata','corr')
else
    %load([bbio_dir,'Stat_',bbio,'_',dataset,'_',kind,'_',parameter,'_',cruise,'_',year,'.mat'])
end

% ------------------------------------------------------
% ---- FIGURES
% ------------------------------------------------------

if (fig==1)
    f_overview(dat,nmodel,dataset,kind,parameter,cruise,bbio,bbio_dir,year)
end


% -----------------------------------------------------------------
% ---- PRIMARY PRODUCTION VALIDATION by BIOGEOCHEMICAL PROVINCES
% -----------------------------------------------------------------

%PP_provinces




